#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
#Copyright 2011. Jason Weirather.
#This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 or any later version.
# two inputs the fasta file since we need all the lengths of the proteins
# and the blast results.

my $input_fasta_name = shift @ARGV;
my $input_results = shift @ARGV;

# if its a protein remember to chop it
# may want the option of showing relation between like groups
# or between out groups.  ie matrix style.

###################################
#       green  red   blue   orange  
#green   png1  png2  png3   png4
#
#red     png5  png6  png7   png8
#
#blue    png9  png10 png11  png12
#
#orange  png13 png14 png15  png16
#
####################################

#those names can change but it will have to be something like that

#first bring in that fasta file 
my $infas = Bio::SeqIO->new('-format' => 'fasta', '-file' => "$input_fasta_name");
my %lengths;
my $is_protein_count = 0;
while( my $seq = $infas->next_seq()) {
  if($seq->seq =~ /[RDEQHILKMFPSWYV]/) { $is_protein_count++; }
  $lengths{$seq->id} = $seq->length;
  my $lastseq =  $seq->subseq($seq->length,$seq->length);
  if($lastseq eq '*') { $lengths{$seq->id}--; } # dont count the star, blast probably doesnt
  #print $lastseq . " $is_protein_count\n";
}
#probably do this simpler.  somehow.
#thats probably good to get the approrpiate lengths.  I hope.
# there may be problems if theres a lot of unknown content.

open(IF,$input_results) or die "no result file.\n";
chomp(my @lines = <IF>);
close IF;

shift @lines;
shift @lines;
shift @lines;
shift @lines;
shift @lines;
#removed the header from the file
#my $types;
#foreach my $line (@lines) { #swing through and find types for each name
#  my @fields = split(/\t/,$line); 
#  $types{$fields[0]} =  $fields[1];
#  $types{$fields[2]} = $fields[3];
#}
my $steps = 20;
my %datasets;
foreach my $line (@lines) {
  my @fields = split(/\t/,$line);
  # each of these is a HSP.  For each of these say whether it is
  # selftype-selftype
  # selftype-targettype
  #
  # also say the length
  # also include a binary map from 1 to 100 showing where within the length
  # of the strand the pairing is occuring
  my %data;
  my @types = sort {$a cmp $b} ($fields[1],$fields[3]); #alphbetize to minimize the number of supertypes 
  my $supertype = "$types[0]---$types[1]";
  #print "$supertype\n";
  #print "$fields[0]\t$fields[2]\t$lengths{$fields[0]}\t$lengths{$fields[2]}\n";
  # figure the histogram for the query and the hit
  my $query_length = $lengths{$fields[0]};
  my $query_type = $fields[1];
  my $hit_length = $lengths{$fields[2]};
  my $hit_type = $fields[3];
  my $query_coverage = $fields[10];
  my $hit_coverage = $fields[12];
  #make an array of zeros for the length of the query and hit
  #make it 20 long since we are only considering histogram distribution of hits
  #that 20 may become an option later.
  my @query_binary_array;
  my @hit_binary_array;
  my $query_step = $query_length / $steps;
  my $hit_step = $hit_length / $steps;
  #now calculated scaled coverage
  my $scaled_query_coverage = scale_coverage($query_coverage,$query_step);
  #print "\t$query_coverage\t$scaled_query_coverage\n";
  my $scaled_hit_coverage = scale_coverage($hit_coverage,$hit_step);  
  my @array_query_coverage = array_of_coverage($scaled_query_coverage,$steps);
  my @array_hit_coverage = array_of_coverage($scaled_hit_coverage,$steps);
  #print "$query_coverage\t$hit_coverage\n";
  #print "$scaled_query_coverage\t$scaled_hit_coverage\n";
  #print "@array_query_coverage\n";
  #print "@array_hit_coverage\n";
  $data{'query_name'} = $fields[0];
  $data{'hit_name'} = $fields[1];
  $data{'query_type'} = $query_type;
  $data{'hit_type'} = $hit_type;
  $data{'query_length'} = $query_length;
  $data{'hit_length'} = $hit_length;
  $data{'query_coverage'} = \@array_query_coverage;
  $data{'hit_coverage'} = \@array_hit_coverage;
  $data{'supertype'} = $supertype;
  # this super type is what it can be compared with
  if(exists($datasets{$supertype})) { push @{$datasets{$supertype}}, \%data; }
  else {
    my @datas;
    push @datas, \%data;
    $datasets{$supertype} = \@datas;
  }
}
#now for each of those data sets you can construct histograms
my @super_histo;
for (my $k = 0; $k < $steps; $k++) { push @super_histo, 0; }
my %histograms;
foreach my $supertype (keys %datasets) {
  #print "$supertype\n";
  my ($type1,$type2) = split(/---/,$supertype);
  if ($type1 eq $type2) { # for when they are the same type make one histogram
    my @histogram;
    for (my $i = 0; $i < $steps; $i++) { push @histogram, 0; }
    my @datas = @{$datasets{$supertype}};
    my $count = 0;
    foreach my $data_ref (@datas) {
      $count = $count + 2;
      my %data = %{$data_ref};
      my @query_coverage = @{$data{'query_coverage'}};
      my @hit_coverage = @{$data{'hit_coverage'}};
      for(my $i = 0; $i < $steps; $i++) { 
        $histogram[$i] = $histogram[$i] + $query_coverage[$i] + $hit_coverage[$i]; 
        $super_histo[$i] = $super_histo[$i] + $query_coverage[$i] + $hit_coverage[$i];
      }
    }
    my %types;
    my %contents;
    $contents{'histograms'} = \%types;
    my %counts;
    $counts{$type1} = $count; #nodes contributing to this distribution (hsp * 2 approx)
    $contents{'node_counts'} = \%counts;
    $types{$type1} = \@histogram;
    $histograms{$supertype} = \%contents;
    #print "@histogram\n";
  } else { #when there are two different types
    my @histogram1;
    my @histogram2;
    my $counts1 = 0;
    my $counts2 = 0;
    for (my $i = 0; $i < $steps; $i++) {
      push @histogram1, 0;
      push @histogram2, 0;
    }
    my @datas = @{$datasets{$supertype}};
    foreach my $data_ref (@datas) {
      my %data = %{$data_ref};
      my @coverage1;
      my @coverage2;
      $counts1++;
      $counts2++;
      if($data{'query_type'} eq $type1) {
        @coverage1 = @{$data{'query_coverage'}};
        @coverage2 = @{$data{'hit_coverage'}};
      } else {
        @coverage1 = @{$data{'hit_coverage'}};
        @coverage2 = @{$data{'query_coverage'}};
      }
      for(my $i = 0; $i < $steps; $i++) { 
        $histogram1[$i] = $histogram1[$i] + $coverage1[$i]; 
        $histogram2[$i] = $histogram2[$i] + $coverage2[$i];
        $super_histo[$i] = $super_histo[$i] + $coverage1[$i]+$coverage2[$i];
      }
    }
    my %types;
    my %contents;
    $contents{'histograms'} = \%types;
    my %counts;
    #print "$type1\t$type2\n";
    $counts{$type1} = $counts1;
    $counts{$type2} = $counts2;
    $contents{'node_counts'} = \%counts;
    $types{$type1} = \@histogram1;
    $types{$type2} = \@histogram2;
    $histograms{$supertype} = \%contents;
    #print "$type1\t@histogram1\n";
    #print "$type2\t@histogram2\n";
  }      
}
#now we have histograms in each of the supertypes, lets normalized them.
my $keynum = scalar(keys %histograms);
open(FO,">$input_results.histo") or die "couldnt write out to .histo\n";
print FO "superkeys: $keynum\n";
foreach my $supertype (keys %histograms) {
  print FO "supertype:$supertype\n";
  my @tys =  keys(%{$histograms{$supertype}->{'histograms'}});
  #print "@tys\n";
  foreach my $type (@tys) {
    print FO "  type:$type\n";
    my @histo = @{$histograms{$supertype}->{'histograms'}->{$type}};
    my $count = $histograms{$supertype}->{'node_counts'}->{$type};
    my @norm_histo = normalize_histo(@histo);
    print FO "  histo:@histo\n";
    print FO "  norm:@norm_histo\n";
    print FO "position\t$type($count)\n";
    for(my $i=0;$i<scalar(@norm_histo);$i++) {
      print FO "$i\t$norm_histo[$i]\n";
    }
  }
}
print FO "combined_histogram\n";
my @norm_super_histo = normalize_histo(@super_histo);
for(my $i=0; $i<scalar(@norm_super_histo); $i++) {
  print FO "$i\t$norm_super_histo[$i]\n";
}
close FO;


sub normalize_histo {
  my @histo = @_;
  #find high and low.
  my $low = 99999999999999;
  my $high = 0;
  foreach my $histo (@histo) {
    if($histo > $high) { $high = $histo; }
    if($histo < $low)  { $low  = $histo; }
  }
  my $range = $high - $low;
  for(my $i = 0; $i < scalar(@histo); $i++) {
    $histo[$i]=($histo[$i]-$low)/$range;
  }
  return @histo;
}

sub array_of_coverage {
  my $coverage = shift @_;
  my $steps = shift @_;
  # make array of zeros for each step
  my @binary;
  for (my $i = 0; $i < $steps; $i++) { push @binary, 0; }
  my @coverages = split(/\s+/,$coverage);
  foreach my $cover (@coverages) {
    my ($start,$finish) = split(/-/,$cover);
    for my $i ($start..$finish) {
      my $index = $i-1;
      $binary[$index] = 1;
    }
  }
  return @binary;
}

sub scale_coverage {
  my $coverage = shift @_;
  my $step_size = shift @_;
  my @coverages = split(/\s+/,$coverage);
  my $n_coverage = "";
  foreach my $cover (@coverages) {
    my ($start,$finish) = split(/-/,$cover);
    #my $scaled_start = $start/$step_size;
    #$scaled_start = sprintf "%.0f", $scaled_start;
    #my $scaled_finish = $finish/$step_size;
    #$scaled_finish = sprintf "%.0f", $scaled_finish;
    my $scaled_start = int($start/$step_size);
    my $scaled_finish = int($finish/$step_size);
    $n_coverage .= "$scaled_start-$scaled_finish ";
  }
  return $n_coverage;
}

=begin MYTEXT
sub getHistoCoverage {
  # oops, we need nucleotide and protein lengths for all these
  # to know how big of a query/hit we are aligning.
  open LEN, 'nt_protein_length.txt';
  chomp(my @lines = <LEN>);
  close LEN;
  my %protein_length;
  my %nucleotide_length;
  foreach my $line (@lines) {
    my @fields = split(/\t/,$line);
    $protein_length{$fields[0]} = $fields[2];
    $nucleotide_length{$fields[0]} = $fields[1];
  }
  my @connections = @_;
  my @hist;
  for (1..100) { push @hist,0; }
  ## decide if its nucleotide or protein based on average last
  ## yeah, i know thats messed up but its the easiest way to code this in right now
  my @lengths;
  my %VSGs;
  foreach my $con (@connections) {
      my $index = $con->{'Index_query'};
      $VSGs{$con->{'node'}} = 1;
      $VSGs{$con->{'target'}} =1;
      my $eee;
      my $bbb;
      foreach my $in (split(/\s/,$index)) {
        ($bbb,$eee) = split(/-/,$in);
      }
      push @lengths, $eee;
  }
  my $stats = Statistics::Descriptive::Full->new();
  $stats->add_data(@lengths);
  my $av_len = $stats->mean();
  my $std_dv = $stats->standard_deviation();
  my $count = $stats->count();
  my $numVSGs = scalar(keys(%VSGs));
  print "avg length is $av_len std_dv is $std_dv for $numVSGs VSGS, $count connect$
  my $isProtein = 1;
  if($av_len > 400) { $isProtein = 0; }
  foreach my $con (@connections) {
   my $index;
   my $unaligned;
   foreach my $which ('query','hit') {
    if($which eq 'query') {
      $index = $con->{'Index_query'};
      $unaligned = $con->{'UnAlign_query'};
    } else {
      $index = $con->{'Index_hit'};
      $unaligned = $con->{'UnAlign_hit'};
    }
    #use both index from the query
    # first is one, last is harder
    my @covers = split(/\s/,$index);
    my $unsofar = 0;
    my $enn = 0;
    my $begg = 0;
    for(my $i = 0; $i < scalar(@covers); $i++) {
      # for each of the gaps
      my $oldenn = $enn;
      ($begg, $enn) = split(/-/,$covers[$i]);
      $unsofar += $begg-$oldenn-1;
    }
    #print "unalign: $unaligned\n";
    #print "unsofar: $unsofar\n";
    #print "enn:     $enn\n";
    #my $last = $unaligned - $unsofar + $enn;
    my $last;
    if($isProtein == 1) {
      if($which eq 'query') { $last = $protein_length{$con->{'node'}}; }
      else { $last = $protein_length{$con->{'target'}}; }
      #print "$con->{'hit'}\n";
      #die;
    } else {
      if($which eq 'query') { $last = $nucleotide_length{$con->{'node'}};}
      else { $last = $nucleotide_length{$con->{'target'}}; }
    }
    #great, now you know the last so you can make array of size 100, and
    #if a particular hit has an alignment on any of the percents you add one to it
    my $interval = $last/100;
    #print "last: $last\n";

    for(my $i = 0; $i < scalar(@covers); $i++) {
      #print "       $covers[$i]\n";
      my ($beg,$en) = split(/-/,$covers[$i]);
      for(my $j = $beg; $j < $en; $j++) { $hist[int($j/$interval)]++;}
    }
   }
  }
  return @hist;
}
=end MYTEXY
=cut


